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Abstract. This paper presents an algorithm for checking temporal prece- 
dence properties of nonlinear switched systems. This class of properties 
subsume bounded safety and capture requirements about visiting a se- 
quence of predicates within given time intervals. The algorithm handles 
nonlinear predicates that arise from dynamics-based predictions used in 
alerting protocols for state-of-the-art transportation systems. It is sound 
and complete for nonlinear switch systems that robustly satisfy the given 
property. The algorithm is implemented in the Compare Execute Check 
Engine (C2E2) using validated simulations. As a case study, a simplified 
model of an alerting system for closely spaced parallel runways is con- 
sidered. The proposed approach is applied to this model to check safety 
properties of the alerting logic for different operating conditions such as 
initial velocities, bank angles, aircraft longitudinal separation, and run- 
way separation. 


1 Introduction 

Dynamic analysis presents a scalable alternative to static analysis for models 
with nonlinear dynamics. The basic procedure for dynamic safety verification 
has three building blocks: (a) a simulation engine, (b) a generalization or bloat- 
ing procedure, and (c) a satisfiability checker. The simulation engine generates 
validated simulations of the model with some rigorous error bounds. The gen- 
eralization procedure uses additional model information to over-approximate 
bounded-time reach sets from the simulations. This additional model informa- 
tion could be, for example, statically computed Lipschitz constants [11], con- 
traction metrics [7] or more general designer-provided annotations [4], Finally, 
the approximation is checked by a satisfiability procedure for inferring safety 
or for iteratively refining its precision. With these three pieces it is possible to 
design sound and relatively complete algorithms for bounded time safety ver- 
ification that also scale to moderately high-dimensional models [4]. 

This paper proposes a new algorithm that extends the reach of the above 
procedure in two significant ways. First, the new algorithm verifies temporal 


2 


precedence properties which generalize bounded safety. A model A satisfies tem- 
poral precedence Pi -<t P 2 if along every trajectory of A, for any time at which 
the predicate If holds, there exists an instant of time, at least b time units sooner, 
where the predicate P\ must hold. The key subroutine in the new verification 
algorithm uses a simulation-based reach set approximation procedure for es- 
timating the time intervals over which the predicates P\ and If may or must 
hold. These estimates are constructed so that the algorithm is sound. The al- 
gorithm is guaranteed to terminate whenever A satisfies the given property 
robustly (relatively complete). That is, not only does every trajectory £ satisfy 
Pi -<b P' 2 , but any small time-shifts and value perturbations of £ also satisfy 
Pi -<(, If. Such relative completeness guarantees usually have the most preci- 
sion that one can hope for in any formal analysis of models involving physical 
quantities. 

Secondly, a new approach to checking satisfiability of nonlinear guarantee 
predicates [8] is proposed. If Pi and P 2 in the above type of temporal prece- 
dence property are in propositional logic or uses linear arithmetic, then exist- 
ing solvers can efficiently check whether a set of states satisfy them. On the 
other hand, if they are written as 3t > 0, f p {x,t ) > 0, where f p is a nonlinear 
real-valued function, then the options are limited. Quantifier elimination is an 
expensive option, but even that is feasible only if f p has a closed form definition 
of a special form (such as polynomial functions). If f p is implicitly defined as the 
solution of a set of ODEs with no analytical solution then quantifier elimination 
is impossible. This paper provides a sound and relatively complete procedure 
for checking bounded time guarantee predicates using simulation-based over- 
approximations of 

The new algorithm is used in the analysis of an interesting and difficult 
verification problem arising from a parallel landing protocol. The Simplified 
Aircraft-based Paired Approach (SAPA) [5] is an advanced operational concept 
that enables dependent approaches in closely spaced parallel runways. In the 
presence of blundering aircraft, the SAPA procedure relies on an alerting algo- 
rithm called Adjacent Landing Alerting System (ALAS) [10]. ALAS uses linear 
and nonlinear projections of the current aircraft positions, velocity vectors, and 
bank angles to detect possible conflicts between the landing aircraft. Given the 
nonlinear characteristics of the ALAS logic, finding operating conditions under 
which the SAPA / ALAS protocol satisfy some safety properties is a challenging 
problem. 

This paper presents a simplified model, written as a switched system, of the 
SAPA / ALAS protocol. The safety properties that are considered on this model 
state that an alert is issued at least b seconds before an unsafe scenario is en- 
countered. These properties are specified as temporal precedence properties of 
the form Alert -<b Unsafe. The proposed verification algorithm is applied to this 
model to formally check these kinds of safety properties for various aircraft and 
runway configurations. 
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2 System Models and Properties 

For a vector v in W, v stands for /: 2 -norm. Given intervals 1,1' , the relation 
I < V holds iff Vi> £ /, W £ V . v < v' . For a realnumber b, I—b = {v — b\v £ I}. 
Subtraction operation over intervals is defined as, I— I 1 = { v — v' | v £ / . v' £ /'}. 
/x /' = {vx |ri £ /, v' £ /'}. For 5 £ R>o and x £ R", B$(x) C R” is the closed 
ball with radius S centered at x. For a set S C R n , Bs(S ) = U x6 si? 5 (r>). For any 
function V : R" x R" £ R> 0 , given a 5 > 0, Bf {x) = {y \ V ( x , y) < 5}. For a set 
S C R 71 , Bj (S') = U xesBf ( x ). For a bounded set A, dia(A) = sup^, |a; — y\ 
denotes the diameter of A. 

A real-valued function a : R>o K>o is called a class K. function if a(0) = 0 
and a is strictly increasing. It is a class /Coo function if additionally a{x) —> oo 
as x — > oo. For a function h : R>o — > R" and a positive real <5 > 0, the (5-left 
shift of h is the function h$ : R>o — > M™ defined as hg(t) = h(t + S) for any 
t £ R>o- A 5-perturbation of h is any function g : R>o —> R n such that for all t, 
| g(t) — h(t) | < 5. A cadlag function is a function which is continuous from the right 
and has a limit from the left for every element in its domain. 

2.1 The Switched System Model 

This paper uses the sivitch system formalism [6] for modeling continuous sys- 
tems. The evolution of an n dimensional switched system is specified by a col- 
lection of ordinary differential equations (ODEs) also called as modes or locations 
indexed by a set I and a szvitching signal that specifies which ODE is active at a 
given point in time. Fixing a switching signal and an initial state, the system is 
deterministic. Its behavior is the continuous, piece- wise differentiable function 
of time obtained by pasting together the solutions of the relevant ODEs. The 
symbol I represents the set of modes and n represents the dimension of the 
system with R" as state space. 

Definition 1. Given the set of modes X and the dimension n, a switched system A 
is specified by the tuple (0, T, E), with 

(i) 0 C R n , a compact set of initial states, 

(ii) T = {/,; : R n —> R ra } ie x, an indexed collection of continuous, locally Lipschitz 
functions, and 

(iii) E, a set ofsivitching signals, where each a £ Eisa cadlag function a : R> 0 —> X. 

The semantics of A is defined in terms of its solutions or trajectories. For a 
given initial state xq £ 0 and a switching signal o £ E, the solution or the 
trajectory of the switched system is a function : ®>o —> R", such that: 

= Xq, and for any t > 0 it satisfies the differential equation: 

i*oA*) = fvwiUAt))- (i) 

When clear from context, the subscripts xq and a are dropped from f Under 
the stated locally Lipschitz assumption of the ff s and the cadlag assumption 
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on a, it is well-known that Equation (1) has a unique solution and that indeed 
the trajectory £ is a well-defined function. 

A bounded time switching signal can be represented as a sequence a = 
mo, mi, . . . , rrik where each rrii is a pair in I x R + , with the two components 
denoted by rrii.mode and rrii.time. The sequence define a(t) = rrii.mode for 
all t £ E}=o mj.time, E' =0 mj.time). A set of switching signals £ is repre- 
sented as a switching interval sequence S = qo, q\ .... <jk, where each q :l is a pair 
with qj.mode £ X and qj. range is an open interval in R> 0 . Given a switch- 
ing interval sequence S, the set sig(S) denotes the set of switching signals 
( 7 = mo, mi, . . . ,mfc, such that rrij.mode = qj.mode and mj.time £ qj. range. 
By abuse of notation, a set of switching signals £ and its finite representation S 
with sig(S) = £ are used intechangeably. The expression width(S) denotes 
the size of the largest interval qi. range. The refinement operation of £, de- 
noted as refine(S), gives a finite set of switching interval sequences S such that 
Us'es sig(S') = sig(S) and for each S' £ S, width(S') < width(S)/2. 

2.2 Temporal Precedence with Guarantee Predicates 

A predicate for the switched system A is a computable function P : R" — y {T, _L} 
that maps each state in R" to either T (true) or _L (false). The predicate is said 
to be satisfied by a state x £ K" if P(x) = T. A guarantee predicate [8] P{x) is 
a predicate of the form 3t > 0, f p (x, t) > 0, where f p : R" x R — > R is called a 
lookahead function. A guarantee predicate holds at a state x if there exists some 
future time t at which f p (x. t) > 0 holds. Using a quantifier elimination pro- 
cedure, a guarantee predicate can be reduced to an ordinary predicate without 
the existential quantifier. However, this is an expensive operation, and more 
importantly, it is only feasible for restricted classes of real-valued lookahead 
functions with explicit closed form definitions. Section 3.1 presents a technique 
to handle guarantee predicates with lookahead functions as solutions to nonlin- 
ear ODE. As seen in Section 4, such lookahead functions are particularly useful 
in designing alerting logics such as ALAS. 

Temporal precedence properties are a class of properties specified by a pair of 
predicates that must hold for any behavior of the system with some minimum 
time gap between them. More precisely, a temporal precedence property & is 
written as f = Pi -</, P 2 , where Pi and P 2 are (possibly guarantee) predicates 
and & is a positive real number. The property f = Pi R& / 2 is satisfied by a 
particular trajectory £ of A iff 

Vt 2 > 0, if P 2 (£(t 2 )) then 3ti, 0 < U < i 2 - 6, Pi(£(fi)). (2) 

In other words, along £, predicate Pl should be should be satisfied at least b 
time units before any instance of P 2 is satisfied. A switched system A satisfies 
<j), if every trajectory of A satisfies 0 . With a collection of precedence properties, 
it is possible to state requirements about ordering of some predicates before 
others. 

The property f is said to be robustly satisfied by a system if 3r > 0, S > 0 
such that all t' < t left shifts and all d -perturbations of all trajectories £ satisfy 
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the property. An execution £ is said to robustly violate a precedence property 
Pi -<i, P 2 if there is a time instant t 2 such that P 2 (f(t 2 )) holds, and for some 
S > 0, all (5-perturbations £' of £ and t\ £ (0, t 2 — b), Pi does not hold in £' at 
time t\ . A system is said to robustly violates <i> = I\ -</, P 2 if some execution £ 
(from an initial state) robustly violates (f>. 

3 Simulation-based Verification of Temporal Precedence 

This section presents an algorithm for verifying temporal precedence proper- 
ties of switched systems and establish its correctness. Similar to the simulation- 
based safety verification algorithm presented in an earlier work [4], this algo- 
rithm has three key features: (a) it uses validated simulations for the dynamics 
in T , (b) it requires model annotations called discrepancy functions for the dy- 
namics in in T . Finally, (c) it requires a procedure for checking satisfiability of 
nonlinear guarantee predicates arising from solutions of differential equations. 

For a given initial state xq and an ODE x = f(x. t) which admits a solution 
£, a fixed time-step numerical integrator produces a sequence of sample points 
e.g., x \ , x 2 , , xi £ R" that approximate the trajectory £ Xo at a sequence of 
time points, say f Xo (h),f Xo (2h), ... ,f Xo (l x h). However, these simulations do 
not provide any rigorous guarantees about the errors incurred during numeri- 
cal approximations. Rigorous error bounds on these simulations, which can be 
made arbitrarily small, are required for performing formal analysis. One such 
notion of a simulation for an ODE is defined as follows. 

Definition 2. Consider an ODE x = f(x, t). Given an initial state, xq, a time bound 
T > 0, error bound e > 0, and a time step r > 0, an (xo, T, e, T)-simulation trace is 
a finite sequence (Ri, [f 0 , ti]), (f? 2 , [ti, t 2 ]), [ti-i, ti)) where each Rj C R", 

and tj £ K> 0 such that Vj, 1 < j < l 

(1) tj - 1 < tj, tj — tj - 1 < r,t 0 = 0, and ti = T, 

(2) Vt S [tj-i,tj],f Xo (t) £ Rj, and 

(3) dia(Rj) < e. 

Numerical ODE solvers such as CAPD 3 and VNODE-LP 4 can be used to gen- 
erate such simulations for arbitrary values of r and e using Taylor Models and 
interval arithmetic. 

Definition 3. A smooth function V : M 2n — ► M> 0 is called a discrepancy function 
for an ODE x = f(x, t), if and only if there are functions a,a£ /Coo and a uniformly 
continuous function ft : M 2 ” xR-> R> 0 with j3(xi,x 2 ,t) — > Oas \xi — x 2 \ 0 such 

that for any pair of states xi,x 2 £ R ra : 


a(\xi — x 2 \) < V(xi,x 2 ) < a(\xi - x 2 \) and 
Vt > 0. V(f Xl (t),f X2 (t)) < P(xi,x 2 ,t), 


(3) 

(4) 


3 http : / / capd .ii.uj.edu . pi /index . php 

4 http : / / www . cas .mcmaster . ca/ -nedialk/vnodelp 
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where £ denotes the solution of the differential equation. A tuple (a, a , /3) satisfying the 
above conditions is called a witness to the discrepancy function. 

The discrepancy function provides an upper bound on the distance between 
two trajectories starting from different initial states x± and x^ ■ This upper bound, 
together with a simulation, is used to compute an overapproximation of the set 
of all reachable states of the system from a neighborhood of the simulation. 

For linear and affine dynamics such discrepancy functions can be computed 
by solving semidefinite programs [4]. In [4], classes of nonlinear ODEs were 
identified for which Lipschitz constants, contraction metrics, and incremental 
Lyapunov functions can be computed. These classes are all special instances 
of Definition 3. For the switched systems A with a set of differential equations 
■F = {fihex, a discrepancy function for each /* (namely, V, and its witness 
(a,;,ab /?*)) is required. Using discrepancy function and validated simulations 
as building blocks, a bounded overapproximation of the reachable set for initial 
set 0, set of switching signals S, and time step r can be defined as follows. 

Definition 4. Given an initial set of states 0, switching interval sequence S, dynam- 
ics T , time step r > 0, and error bound e > 0, a (0, S, e, r)-ReachTube is a sequence 
if = (Or, [t 0 ,ti]), (0 2 , [ti,t2]), • • • , (Oi, where Oj is a set of pairs ( R,h ) 

such that R C R", and h £ X, such that, Vj, 1 < j < l 

(1) tj — i ^ tj, tj tj — i C r, to — 0/ 

(2) Vx 0 G <9,Ver G sig(S)fyt G [tj-i,tj],3(R,h) G Oj, such that, £ XOjCr (t) G 
R, a(t) = h, 

(3) V(R,h) G Oj,dia(R) < e, and 

(4) each mode in X occurs at most once in Oj. 

Intuitively, for every given time interval [tj-i,tj], the set O t contains an 
(R, h) pair such that It overapproximates the reachable set for the mode h in 
the given interval duration. In a previous work on verification using simula- 
tions [4], an algorithm that computes overapproximation of the reachable set 
via sampled executions and annotations is presented. The procedure, called 
ComputeReachTube, takes as input the initial set 0, switching signals S, par- 
titioning parameter S, simulation error e' , and time step r. It compute the se- 
quence if and error e such that h is a (0, S, e, r)-ReachTube. The procedure is 
outline below. 

1. Assign to Q, the set of initial states 0. 

2. For each ^ in the switching interval sequence S = r/ 0 . q x . ... . q k . 

3. Compute X = {x\,X 2 , ■ ■ ■ ,x m }, a (5-partitioning of <5, such that <5 C U B$(xi). 

4. Generate a validated simulation (Definition 2) ij for every state x G X with 
error e' , time step r. Then, compute the ReachTube for Bs(x o) by bloating 

p as B^ qi ' mode (77), where e = sup{f qi . rnode (y,x,t)\y G B s (x)}. 

5. Compute the union of each of the ReachTubes for B$(xf) as the ReachTube 
for mode qi.mode. 

6. Compute the initial set for the next mode by taking the projection of ReachT ube 
for qi.mode over the interval qt.range as Q. Repeat steps 3-6 for q l+ \ . 
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The order of overapproximation of the ReachTube computed using the proce- 
dure described above is the maximum bloating performed using the annotation 
Vqi.mode and qi .mode for all the modes in S. This overapproximation and the 
error in simulation gives the value of e such that ^ is a (0, S, e, r)-ReachTube. 
The nondeterminism during the switching times from one mode to another en- 
ables the reachable set to be in two different modes at a given instance of time, 
which is reflected in Oj. The following proposition holds [4]. 

Proposition 1. Given an initial set 0, switching signals S, partitioning parameter S, 
simulation error e' and time step r, let {ip, e) = ComputeReachTube(0, S, 5, e', r). As 
dia{0) —> 0, width{S) — > 0, S — > 0, e' — > 0, and r — > 0, then e — > 0. 

3.1 Temporal Precedence Verification Algorithm 

Check Refine (see Figure 1) performs the following steps iteratively: (1) Create 
an initial partition of the set of start states 0. (2) Compute the ReachTubes for 
each these partitions as given in Definition 4. (3) Check the temporal precedence 
property for the ReachTube. (4) Refine the partitioning if the above check is 
inconclusive, and repeat steps (2)-(4). 

A key step in the procedure is to verify whether a given ReachT ube satisfies 
a temporal precedence property The collection of intervals rrmstlnt, notint, and 
maylnt, for a given predicate P and ReachTube ip are used in the verification of 
temporal precedence properties. They are defined as follows. 

Definition 5. Given a ReachTube ip = (O i, [t 0 , ti]), • ■ ■ , (Oz, t;]) and a pred- 
icate P,for all j > 0, 

[tj_i,tj] £ mustInt{P , ip) iff V(f?, h) £ Oj,RQ P. 

[tj-i,tj\ £ notInt{P,ip) iff W{R,h) £ Oj,R C P c . 

[tj-i,tj] £ mayInt{P,ip) otherwise. 

Definition 5 classifies an interval [tj-i , t.j] as an element of mustInt{P,ip ) 
only if the overapproximation of the reachable set for that interval is contained 
in P. Similar is the case with notint {P, ip). However if the overapproximation 
of the reachable set cannot conclude either of the cases, then the interval is clas- 
sified as mayInt{P,ip). There are two possible reasons for this: first, the order 
of overapproximation is too coarse to prove containment in either P or P c ; sec- 
ond, the execution moves from the states satisfying P to states not satisfying P 
during that interval. Thus, better estimates of mustlnt, notint and maylnt can 
be obtained by improving the accuracy of ReachTube ip. 

To compute mustlnt, maylnt, and notint as defined in Definition 5, it is nec- 
essary to check if R C P or R C P c . However, for guarantee predicates with 
lookahead functions that use the solutions of ODEs, it is unclear how to per- 
form these checks. Section 3.2 describes a simulation-based method to address 
this challenge. The algorithm in Section 3.2 will, in fact, provide weaker guar- 
antees. Assuming P is an open set, the algorithm will answer correctly when 
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R C P and when for some 5 > 0, Bs (R) C P c ; in other cases, the algorithm may 
not terminate. Such weaker guarantees will turn out to be sufficient for the case 
study considered in this paper. 

Definition 6. Given ReachTube ip and temporal precedence property Pi ^ P 2 , ip is 
said to satisfy the property iff for any interval I', I' £ mustInt(P 2 , ip)UmayInt(P 2 , ip), 
exists interval I, I £ mustInt(Pi , ip) such that I < I' — b. Also, ip is said to violate the 
property if 31' £ mustInt(P 2 ,ip) such that, V/ £ mustInt(Pi,ip) U mayInt(Pi,ip), 

V -b< I. 

From Definition 6 it is clear that if a ReachT ube ip satisfies a temporal prece- 
dence property, then for all the trajectories corresponding to the ReachTube, the 
predicate Pi is satisfied at least b time units before P 2 . Also, if the ReachTube 
violates the property, then it is clear that there exists at least one trajectory such 
that for an instance of time, i.e., in /' £ mustInt(P 2 ,ip) at all the time instances 
at least b units before, the predicate P\ is not satisfied. In all other cases, the 
ReachTube cannot infer whether the property is satisfied or violated. As this 
inference depends on the accuracy of mustlnt, notint and maylnt. More accu- 
rate ReachT ubes produce better estimates of these intervals and hence help in 
better inference of temporal precedence property. 

Given a system A and property Pi -p, P 2 , one can compute the ReachTube 
for the system and apply Definition 6 to check whether the system satisfies 
the temporal precedence property. This is however not guaranteed to be useful 
as the approximation of ReachTube computed might be too coarse. The algo- 
rithm CheckRefine refines, at each iteration, the inputs to compute more precise 
ReachT ubes. Proposition 1 guarantees that these ReachT ubes can be made ar- 
bitrarily precise. 

The algorithm (in Figure 1) first partitions the initial set into ^-neighborhoods 
(line 4) and compute ReachT ubes for every switching interval sequence in fl 
(line 7). If all these ReachT ubes (that is all the executions from neighborhood) 
satisfy the property, then the neighborhood is removed from Q. Similarly, the al- 
gorithm CheckRefine returns that the property is violated only when ReachTube 
violates the property. If neither can be inferred, then the parameters to func- 
tion ComputeReachTube are refined in line 11 to increase their precision. Since 
this operation is iteratively performed to obtain arbitrarily precise ReachT ubes, 
Soundness and Relative completeness follow from Definition 6 and Proposi- 
tion 1. 

Theorem 1 (Soundness and Relative Completeness). Algorithm CheckRefine is 

sound, i.e., if it returns that the system satisfies the property, then the property is indeed 
satisfied. If it returns that the property is violated, then the property is indeed violated 
by the system. Further, if predicates Pi and P 2 are open sets, and there is a procedure 
that correctly determines if for a set R, R C P t (for i = 1,2 ) or if there is S > 0 such 
that Bg(R) C Pp (for i = 1,2). Then, if the system A satisfies Pi -<; h P 2 or robustly 
violates Pi -<b P 2 then CheckRefine terminates with the right answer. 
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1: Input: A = (0,J°, E), {V, (a i: a it Pi)} ten, Pi -<b P 2 , S 0/ So, e' 0 , r 0 . 

2: Q 4 — 0; 17 <— {E}; S <— <5o; S' «— So; t <— e' 0 ; r «— to 

3: while Q ^ 0 do 
4: A’ 8-partition(Q); 

5: for all xo € X do 

6: for all S € 17 do 

7: {ip , e) = ComputeReachTube(_Ba(a;o), S', <T, t, r) 

8: if ip satisfies Pi Xt, P 2 then continue; 

9: else if ip falsifies Pi ~ii, P 2 return "Property Pi -<i, P 2 is violated" 

10: else 

11: Q <r- n \ {S'} U refine(S); S <r- 8/2; S' «- 5'/2, e' <- e'/ 2; r <- r/2; 

12: goto Line 4 

13: end if 

14: end for 

15: Q Q \ B s (xo) 

16: end for 

17: end while 

18: return "Property Pi Xb P 2 is satisfied". 

Fig. 1. Algorithm CheckRefine: Partitioning and refinement algorithm for verification of 
temporal precedence properties. 


3.2 Verification of Guarantee Predicates 

As discussed in the Section 2.2, guarantee predicates are of the form P(x ) = 
3t > 0, f p (x , t) > 0, where f p is called a lookahead function. Section 3.1 presents 
an algorithm for time bounded verification of such predicates of the special 
form P{x) = 30 < t < Ti,w p (£' x (t)) > 0, where w p is a continuous function 
and £' is solution of ODE y = g(y. t). The algorihm CheckGuarantee in Fig- 
ure 2 checks whether R C P or an open cover of R is contained in P c has 
been defined. This algorithm, similar to CheckRefine, computes successively 
better approximations for the ReachTube and checks whether the predicate 
P' = w p {x) > 0 is satisfied by the reach tube. This is done by calculating 
mustInt{P' ,ip) and mayInt(P' ,ip) as defined in Definition 5. If the mustlnt is 
non-empty, then it implies that the predicate P is satisfied by the ReachTube 
and hence R C P. If both the maylnt and mustlnt are empty sets, then, clearly 
the predicate P is not satisfied in the bounded time T; by any state in R, and 
hence an open cover of R is contained in P c . Soundness and Relative Complete- 
ness of CheckGuarantee follow from CheckRefine (proofs in full version 5 ). 

4 Case Study: A Parallel Landing Protocol 

The Simplified Aircraft-based Paired Approach (SAPA) is an advanced opera- 
tional concept proposed by the US Federal Aviation Administration (FAA) [5]. 

5 https : / / wiki .cites . illinois . edu/wiki/ di spl ay /Mi traRe search/ 
Verif ication+of +a+Parallel+Landing+Protocol 
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1: Input: R, y = g(y,t), S', V g {xi,x 2 ), (a g ,a g ,/3 g ) w p , 5, r, Ti 

2: while R ^ 0 do 
3: X «— S-partition(R); 

4: for all xo € X do 

5: {ip, e) = ComputeReachTube(_B{ (xo), <S ,/ , <5, S, r); 

6: if mustlnt(w p , ip) ^ 0 then R 4— R\ Bs(x o) 

7: else if mustlnt(w p ,ip) U maylnt(w p , ip) = 0 then return "UNSAT" 

8: end if 

9: end for 

10: 5 <— 5/2; r 4— r /2; 

11: end while 
12: return "SAT". 

Fig. 2. Algorithm CheckGuarantee: Decides whether a lookahead predicate is satisfied in 
a given set ft 


The SAPA concept supports dependent, low-visibility parallel approach opera- 
tions to runways with lateral spacing closer than 2500 ft. A Monte-Carlo study 
conducted by NASA has concluded that the basic SAPA concept is technically 
and operationally feasible [5]. SAPA relies on an alerting mechanism to avoid 
aircraft blunders, i.e., airspace situations where an aircraft threats to cross the 
path of another landing aircraft. 

NASA's Adjacent Landing Alerting System (ALAS) is an alerting algorithm 
for the SAPA concept [10]. ALAS is a pair-wise algorithm, where the two air- 
craft are referred to as oivnship and intruder. When the ALAS algorithm is de- 
ployed in an aircraft following the SAPA procedure, the aircraft considers itself 
to be the ownship, while any other aircraft is considered to be an intruder. The 
alerting logic of the ALAS algorithm consists of several checks including con- 
formance of the onwship to its nominal landing trajectory, aircraft separation at 
current time, and projected aircraft separation for different trajectories. 

A formal static analysis of the ALAS algorithm is challenging due to the 
complexity of the SAPA protocol and the large set of configurable parameters 
of the ALAS algorithm that enable different alerting thresholds, aircraft perfor- 
mances, and runway geometries. This paper considers the component of the 
ALAS alerting logic that checks violations of predefined separation minima 
for linear and curved projected trajectories of the current aircraft states. This 
component is one of the most challenging to analyze since it involves nonlin- 
ear dynamics. Safety considerations regarding communication errors, pilot and 
communication delays, surveillance uncertainty, and feasibility of resolution 
maneuvers are not modeled in this paper. 

For the analysis of the landing protocol, this paper considers a blundering 
scenario where the intruder aircraft turns towards the ownship during the land- 
ing approach. The dynamics of the aircraft are modeled as a switched system 
with continuous variables sxi, syt, vxi, vyi and sx 0 , sy 0 , vx or and vy 0 repre- 
senting the position and velocity of intruder and ownship respectively. The 
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switching system has two modes: approach and turn. The mode approach rep- 
resents the phase when both aircraft are heading towards the runway with con- 
stant speed. The mode turn represents the blundering trajectory of intruder. In 
this mode, the intruder banks at an angle <f>i to turn away from the runway to- 
wards the ownship. The switching signal determines the time of transition from 
approach to turn. In this mode, the differential equation of the ownship remains 
the same as that of approach, but the intruder's turning motion with banking 
angle <j>i is 
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where c x and c y are constant functions of the initial states of the ownship and 
intruder, and w, is the angular speed of intruder. Given the bank angle b,, the 
angular speed is given by Wi = f '? ' a " f ' ' =, where G is the gravitational con- 

° 1 ° j yj VXi 2 -\-Vyi 2 ° 

stant. The upper bound on the bank angle <f>i is denoted as </> max . 

The system starts in the approach mode with the initial position of the in- 
truder at sXi = syi = 0 and the ownship at sx a = xsep and sy Q = ysep, where 
xsep denotes the lateral separation between the runways and ysep denotes the 
initial longitudinal separation between the aircraft. The initial velocities of both 
aircraft along the x-axis are 0 and the initial velocities along the y-axis are pa- 
rameters. The time of switching from approach mode to turn mode is nondeter- 
ministically chosen from the interval T SW i tc h = [2.3, 2.8]. These parameters and 
the initial values of the variables are constrained by the SAPA procedure [5]. 

4.1 Alerting Logic and Verification of Temporal Precedence Property 

The alerting logic of ALAS considered in this paper issues an alert when the air- 
craft are predicted to violate some distance thresholds called Front and Back [10] 
To predict this violation, the aircraft projects the current state of the system with 
three different dynamics: first, the intruder does not turn, i.e., banking angle 0°, 
second, the intruder turns with the specified bank angle <j>i and third, the in- 
truder turns with the maximum bank angle 4> m ax • If any of these projections 
violates the distance thresholds, then an alert is issued. The alert predicates 
for the each one of these projections are represented by Alerts, Alert and 
Alert^max' respectively. Thus, the alerting logic considered in this paper is de- 
fined as Alert = Alertg V Alert ^ V Alert j, max . 

The alert predicates Alert 0 , Alert q>i and Alert,p mar are guarantee predicates. 
The lookahead function for Alerts is defined as follows: from a given state x, 
it computes the projected trajectory of the aircraft when intruder turns at bank 
angle n. If these trajectories intersect, then it computes the times of intersection. 
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That is, it computes C,t 0 such that sir '(C) = sx' 0 (t 0 ) and sy'(C) = sy' 0 {t 0 ), where 
sx ( . sy[, sx' n , sy' 0 represent the positions of the intruder and ownship aircraft in 
the projected trajectory. If such C and t 0 exist, the Alerts is defined as: 

Alert ^(x) = iff C > to ? (At 2 x (vx 2 + vy 2 ) < Back 2 ) 

: (At 2 x (vx 2 + vy 2 ) < Front 2 ), 

where At = tj — t a . If such tj and t D do not exist, then Alert v ( x) = _L. The 
expression al b : c is a short hand for if (a) then b else c. 

As the guarantee predicates cannot be handled by SMT solvers. Section 3.2 
proposes a simulation based algorithm for handling them. In this case study, the 
proposed technique is used to resolve the nonlinearities of t a and t, in the Alerts 
predicate. An overapproximation Alert of Alerts is computed as: Alert' n (x) = 
T if and only if 

Ti > T 0 ? (AT 2 x (vx 2 0 + vy 2 0 ) < Back 2 ) 

: (AT 2 x (vx 2 + vy 2 ) < Front 2 ) 

where AT = Ti —T 0 . The numerical values of Ti and T a computed simplify the 
Alert^ predicate and can be handled by SMT solvers. 

A state of the system where the intruder aircraft is inside a safety area sur- 
rounding the ownship is said to be unsafe. This paper considers a safety area of 
rectangular shape that is SafeHoriz wide, starts a distance SafeBack behinds the 
ownship and finishes a distance SafeFront in front of the ownship. The values 
SafeHoriz, SafeBack and SafeFront are given constants. Formally, the predicate 
Unsafe is defined as Unsafe(x) = (syi > syf^syi — sy 0 < SafeFront : sy 0 — syi < 
SafeBack) and | sxi — sx 0 \ < SafeHoriz. 

The correctness property considered in this paper is that an alert is raised 
at least b seconds before the intruder violates the safety buffer where b is in 
the range [4, 15]. This can written as a temporal precedence property Alert -<i, 
Unsafe. 

4.2 Verification Scenarios and C2E2 Performance 

The verification algorithms of Section 3 are implemented in the tool Compute 
Execute Check Engine (C2E2). C2E2 accepts Stateflow (SF) charts as inputs, 
translates them to C++ using CAPD for generating rigorous simulations. For 
checking SAT queries, it uses Z3 [1] and GLPK 6 . The discrepancy functions for 
the aircraft dynamics were obtained by computing incremental Lyapunov-like 
function using MATLAB [4]. The following experiments were performed on 
Intel Quad Core machine 2.33 GHz with 4GM memory. 

The temporal precedence property Alert -</, Unsafe is checked for several 
configurations of the system, i.e., values of parameters and initial values of state 
variables. For these experiments, the time bound for verification is set to 15 
seconds and the time bound for projection is set to 25 seconds. 

6 http : //www . gnu . org/ software/ glpk 
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(a) Scenario 1. (b) Scenario 2. (c) Scenario 3. 

Fig. 3. Figure depicting the set of reachable states of the system. Color coding is used to 
depict whether the alert is issued by the alerting algorithm 


Scenario 1. The system configuration is specified by the following param- 
eters and variables: xsep £ [0.22,0.24] km, ysep £ [0.2,04] km, fa = 30°, 
4>max = 45°, vy 0 = 0.07 km/s and vyi = 0.08 km/s. With this configura- 
tion, C2E2 proves that the system satisfies the temporal precedence property 
Alert -<4 Unsafe and an alert is generated 4.38 seconds before the safety is 
violated. The set of reachable states of the ownship and the intruder when the 
safety property is violated is shown in red and the safe states reached are shown 
in blue and green respectively in Figure 3(a). 

Scenario 2. Increasing the intruder 
velocity to vyi = 0.11 km/s, and bank 
angle fa = 45° from the configuration 
of Scenario 1 results in Scenario 2. In 
this case, the safe separation between 
the intruder and the ownship is always 
maintained as the intruder completes 
the turn behind the ownship. Also, the 
alarm is not raised and hence the prop- 
erty Alert -<b Unsafe is satisfied. 

Scenario 3. Changing the configu- 
ration by vyi — 0.11 km/s, xsep £ 

[1.02, 1.04] km, and fa = 45° from 

. - . _ . _ lable 1. Running times. Columns 2-5: Verification 

bcenario 1 results in Scenario 3. (_2.hz Result, Running time, # of refinements, value of b for 

proves that the simplified alerting logic which A ~ <h u 15 satisfied. 

considered in this paper issues a false- 

alert, i.e., an alert is issued even when 

the safety of the system is maintained. Though the property Alert -<j, Unsafe is 
not violated, avoiding such circumstances improves the efficiency of the proto- 
col and C2E2 can help identify such configurations. 

Scenario 4. Placing the intruder in front of ownship, i.e., ysep = —0.3 km and 
vyi = 0.115 km/s from configuration in Scenario 1 results in Scenario 4. C2E2 
proves that the simplified alerting logic considered in this paper misses an alert, 
i.e., does not issue an alert before the safety separation is violated. Such see- 
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narios should always be avoided as they might lead to catastrophic situations. 
This demonstrates that C2E2 can aid in identifying scenarios which should be 
avoided and help design the safe operational conditions for the protocol. 

Scenario 5. Reducing the xsep £ [0.15,0.17] km and ysep £ [0.19,0.21] km 
from configuration in Scenario 1 gives Scenario 5. For this scenario, C2E2 did 
not terminate in 30 mins. Since the verification algorithm presented in Sec- 
tion 3 is sound and relatively complete only if the system robustly satisfies 
the property, it is conjectured that Scenario 5 does not satisfy the property ro- 
bustly. Upon closer inspection, it is observed that the partitioning parameter 
S = 0.0005 and time step r = 0.001 (typical values at termination are S = 0.005 
and r = 0.01), which support the conjecture. 

The running time of verification procedure and their outcomes for several 
other scenarios are presented in Table 1 . Scenarios 6-8 introduce uncertainty in 
the initial velocities of the aircraft with all other parameters remaining the same 
as in Scenario 1. The velocity of the aircraft are changed to be vy 0 £ [0.07, 0.075] 
in scenario 5, vyi £ [0.107, 0.117] in scenario 6, and vXi £ [0.0, 0.005] in scenario 
7 respectively. Scenarios x.l is similar to Scenario x (for x being 6,7,8), however 
with twice the uncertainty in the velocity. Scenario 9 is obtained by changing 
the runway separation to be xsep = 0.5 ±0.01. Scenario 10 is obtained by reduc- 
ing the xsep = 0.2 ± 0.01. Scenario x.l is similar to Scenario x (for x being 9,10) 
however with twice the time horizon for verification and projection. These re- 
sults suggest that the verification time depends on time horizon approximately 
linearly. 


5 Related Work and Conclusion 

There are several MATLAB based tools for analyzing properties of switched 
systems using simulations. Breach [3] uses sensitivity analysis for analyzing 
STL properties of systems using simulations. This analysis is sound and rela- 
tively complete for linear systems, but does not provide formal guarantees for 
nonlinear systems. S-Taliro [9] is a falsification engine that search for counterex- 
amples using Monte-Carlo techniques and hence provides only probabilistic 
guarantees. STRONG [2] uses robustness analysis for coverage of all execu- 
tions from a given initial set by constructing bisimulation functions. Currently 
this tool computes bisimulation functions for only linear or affine hybrid sys- 
tems and does not handle nonlinear systems. 

This paper presents a dynamic analysis technique that verifies temporal prece- 
dence properties and an approach to verify guarantee predicates that use solu- 
tions of ODEs as lookahead functions. These techniques are proved to be sound 
and relative complete. The verification approach is applied to a landing proto- 
col that involves nonlinear dynamics. The case study demonstrated that the 
proposed technique can not only verify safety properties of the alerting logic, 
but also could identify conditions for false and missed alert which are crucial 
in designing the operational concept. 
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